######################################################
# CBQ functionality
######################################################

library(rstan)
rstan_options(auto_write = TRUE)
options(mc.cores = parallel::detectCores())

pald = function(x,p){
  out = ifelse(x<0,p*exp(x*(1-p)), 1-(1-p)*exp(-x*(p)))
  return(out)
}

cbq <- function(N, # number of obs
               n_covariate, # number of covariates
               X, # covariates
               Y, # d.v
               qtl, # quantile to be estimated 1...9
               nchain = 2, # number of parallel chains
               niter = 1000, # number of interations
               thin = 1, # thinning
               seed = 1234567, # random seeds
               offset = 1e-20, # offset to enhance numerical stability
               max_treedepth = 10
               ){
  datlist = list(N=N,
                 Y = Y,
                 D = n_covariate,
                 X = X)
  pars = c("alpha","beta")
  ms_cbq_model = stan_model(paste("replication_Rasmussen/stan_new_models/cbq_binary_q",qtl,".stan",sep=""))
  ms_cbq_stan = sampling(ms_cbq_model, 
                         data=datlist,
                         pars=pars, 
                         seed=seed,
                         chains=nchain,
                         iter=niter,
                         thin=thin,
                         control = list(max_treedepth = max_treedepth))
  return(ms_cbq_stan)
}
